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ABSTRACT 

We have studied formation of planetesimals at a radial pressure bump in 
a protoplanetary disk created by radially inhomogeneous magnetorotational in- 
stability (MRI), through three-dimensional resistive MHD simulations including 
dust particles. In our previous papers, we showed that the inhomogeneous MRI 
developing in non-uniform structure of magnetic field or magnetic resistivity can 
transform the local gas flow in the disk to a quasi-steady state with local rigid 
rotation that is no more unstable against the MRI. Since the outer part of the 
rigid rotation is super-Keplerian flow, a quasi-static pressure bump is created 
and dust concentration is expected there. In this paper, we perform simulations 
of the same systems, adding dust particles that suffer gas drag and modulate gas 
flow via the back-reaction of the gas drag (dust drag). We use ~ 0(10'') super- 
particles, each of which represents ~ O(10^)-O(10^) dust particles with sizes of 
centimeter to meter. We have found that the dust drag suppresses turbulent 
motion to decrease the velocity dispersion of the dust particles while it broadens 
the dust concentrated regions to limit peaky dust concentration, compared with 
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the simulation without the dust drag. We found that the positive effect for the 
gravitational instability (reduction in the velocity dispersion) is dominated over 
the negative one (suppression in particle concentration). For meter-size particles 
with the friction time Tf ~ where VL is Keplerian frequency, the gravita- 
tional instability of the dust particles that may lead to planetesimal formation is 
expected. For such a situation, we further introduced the self-gravity of dust par- 
ticles to the simulation to demonstrate that several gravitationally bound clumps 
are actually formed. Through analytical arguments, we found that the planetes- 
imal formation from meter-sized dust particles can be possible at ~ 5AU, if dust 
spatial density is a few times larger than that in the minimum mass solar nebula. 

Subject headings: protoplanetary disks — instabilities — MHD — planetary 
systems: formation — turbulence 



1. Introduction 



Planets form from coalescence of planetesimals in a protoplanetary disk. Planetesimals 
with more than kilometer sizes should form from dust grains that are initially less than 
micrometer sizes. However, so called "meter-size barrier" exists. Because meter-size particles 
are marginally coupled with disk gas motion and the disk gas rotates slightly slower than 
Keplerian motion due to radially negative pressure gradient of the gas, the particles suffer 
"headwind" and rapidly migrate towa rd the host star . The infall timescale is onl y a few 
hundred years for meter-size particles (jWeidenschillingj 119771 : iNakagawa et al.l Il98ll ) , which 
is much shorter than the growth timescale of particles by mutual collisions. It has not been 
clarified how the particles grow over meter-sizes before infalling to the host star. 

One way to bypass the meter-size barrier is to form clumps f rom dust particles through 

self-g r avit at ional instability (Gl), which occurs on orbital periods ( Safronovlll969 : Goldreich fc Ward 
19731 ). if the dust particles locally have a large enough spatial density. Once bodies of 
kilometer-size or more are formed, they no longer undergo rapid migration. Original idea 
for dust concentration for onset of Gl was vertical settling of dust grains onto the disk mid- 
plane. However, the dust settling induces Kelvin-Helmholtz instability due to difference in 
rotation velocities between the dust-rich layer (Keplerian) and an overlaying dust-poor layer 
(sub-Keplerian), and it prevents the dust layer from becoming dense enough for Gl (e.g 



Weidenschillindll980l : ICuzzi et al.lll993l : ISekivalll998l : llshitsu fc Sekivall2003l: 



unless initial dust to gas ratio in the disk is sufficiently high (e.g., IChiang 
2010al lb[). 



BarrancolEoogj) 



20081: 



Lee et al. 
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Other than the induced KH instabihty, global turbulence is likely to exist in the disk. 
While the turbulence gener ally scatters dust particles, it could concentrate dust particles in 



anti-cyclonic vortexe s (e.g.. iBarge &: Sommerialll995l : IChavanis II2OOOI : iJohansen et al.l 12004 



Inaba fc Barg^boooh . For this mechanism to lead to GI, relatively high initial surface den- 
sity of dust and very weak turbulence may be required. In the case of strong turbulence, 
dust particles have too high velocity dispersion for GI and th e high collision velo c ity between 



dust p articles results i n fragm entation rather than growth ( iGiittler et al.ll2009t IZsom et al. 



2OIOI ). iJohansen et al.l (120071 ) performed local three-dimensional MHD simulation including 
dust particles and showed that weakly fluctuating pressure bumps are created by magne- 
torotational instability (MRI) and meter-size bodies are concentrated at the bumps. In the 
relatively-weak turbulence with the viscosity a ~ 10^^, the dust particles could stay long 
enough and increase their density to cause GI. They found that back-reaction of drag force 
from gas to the dust particles, which we hereafter call "dust drag force," modulates gas 
motion to follow the particles in dust-accumulated regions and weaken the turbulence. 

Although vertical sedimentation of dust is inhibited by KH instability, radial accumula- 
tion is possible. For example, radial dependence of speed of dust migr ation due to gas drag 



can enhance the dust to gas ratio to facilitate GI in inner disk regions (lYoudin fc Shul 12002 



Youdin fc Chiangl 120041 ). This radial migration induces "streaming instability" if dust drag 
is considered. In local dust-accumulated areas (dust clumps), the dust drag force modulates 
the gas flow closer to Keplerian rotation. As a result, "head wind" to the clumps becomes 
weaker and their radial migration due to the gas drag becomes slower, which leads to rapid 



growth of the clum ps by capturing dust partic 
outer regions (e.g., 



Youdin fc Goodman 12005 



2007: Johansen et al. 



2009 



1^ 



es and smaller clumps that migrate faster from 



Youdin fc Johansen 2007; Johansen fc Youdin 



Bai fc StonelboiOaU blJch. The suppression of local turbulence by 



the dust drag decreases the velocity dispersion of the dust particles in the clumps, which is 
also favorable for the GI. 

A global radial pressure bump also leads to radial concentration of dust. The inner 
boundary of "dead zone" is one of such locations. The growt h rate of MRI depends on 



the m agnetic strength and the ionization degree of disk gas (e.g.. lJirull996l : ISano &: Miyama 
I999I ). In the region where the gas ionization degree is low enough or the vertical magnetic 
field is weak enough, the ohmic dissipation decays MRI there ("dead zone"). The disk gas 



is ionized by thermal ioniza tion, the stellar X-rays fe.g . Ilgea &: Glasseo. 



rays (e.g.. 



Umebavashilll983f) and the radionuclides (e.g. IStepinskilll992l ). 



d 1999), the cosmic 



Gammid ( 119961 ) and 



Sano et al.l (120001 ) showed that the dead zone exists in the disk and it is confined in the inner 
disk (< lOAU) near the disk midplane. Because the viscosity is lower in the dead zone and 
disk accretion flux is conserved between the dead and active zones, disk gas column density 
is enhanced in the dead zone. The positive radial gradient of gas column density at the inner 



-4- 



boundary of the dead zo ne produces a pressure bump in which dust particles accumulate 
( iDzyurkevich et al.ll2010l ). However, the inner boundary may be located at < lAU, so the 
planetesimal formation there may be unable to account for formation of icy planets and cores 
of gas giants. 

The column density of tiny dust grains may be enhanced around a snow line due to slow- 
down of the dust radial migration speed by diffusion of sublimated vapor (ICuzzi fc Zahnle 
20081 ) or by down-sizing throu gh sublimation of a icy mantle around a silicate core of dust 
particles ( Saito fc Sironoll2011 ). Since the ionization degree depends on the abundance of the 



tiny grains ( Sano et al. 200ol ). the ioni zation degree would significantly dec rease around the 



Kretke & Lin 


2007; 


Ida & Lin 


2008) 



of the local dead zo n e is a favorable sit e for rapid dust growth in relatively o uter regions 



Brauer et al.ll2008l ). iKato et al.l ( 120091 . which are referred to as Paper I) and iKato et al. 



2010l . Paper II) pointed out that if the local dead zone induced by the snow line is embed- 



ded in the global MRI active zone, the divided inner active zone is sandwiched by the inner 
global dead zone and the outer local dead zone and it can be a stable barrier for dust radial 
migration in which dust particles are accumulated. Even if the snow line is not in the global 
active zone, near the outer boundary of dead zone where MRI is marginal, fluctuations of 
magnetic fields and/or ionization degree could make radially nonuniform MRI structure that 
can be approximated by an active zone radially sandwiched by dead zones. 

In Paper I and II, we have investigated evolution of gas flow of the active zone radially 
sandwiched by dead zones, through shearing-box magnetohydrodynamic simulations. In 
Paper I, performing two-dimensional simulations, we found that the angular velocity profile 
of gas is modified by local MRI turbulence in radially non-uniform magnetic field. The 
vigorous angular momentum and mass transport associated with the MRI turbulence lead 
to a local rigid rotation in the originally active zone. The MRI turbulence can decay to the 
viscosity level of a ~ 10~^ after the transformation to the quasi-steady state with the local 
rigid rotation, because there is no shear motion to create MRI while magnetic field remains. 
In the outer part of the local rigid rotation, gas rotation is super-Keplerian and a pressure 
bump is formed. Note that this gas flow structure is stable. If the rigid rotation is broken, 
the induced shear motion again produces MRI and transports angular momentum and mass 
to recover the rigid rotation as long as the strong enough magnetic field remains. 

In Paper II, we found the same local rigid rotation in the three-dimensional (instead 
of two-dimensional) simulation as shown in Figure [T]d. The calculations with test particles 
show that the boundary region between sub- and super-Keplerian zones acts as a strong and 
stable barrier for the dust migration and it leads to dust particle concentration up to 10,000 
times of the initial value (Figured]:), which could eventually lead to planetesimal formation. 
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However, the dust drag force onto gas was neglected in Paper II, although it would affect 
the dust concentration and velocity dispersion of dust particles in the concentrated regions 
by altering the gas flow. The drag lowers velocity dispersion of dust particles to facilitate the 
GI, while it broadens the dust accumulated region and suppresses peaky dust concentration 
that is rather negative for the GI. The latter effect is positive for the GI, while the former is 
negative. In this paper, we include the dust drag force to the simulations. The equations and 
initial setup employed in our simulation are described in section |2l In section [3], we show the 
simulation results. In section HJ we estimate the possibility of the planetesimal formation 
by analytical arguments. We also demonstrate the planetesimal formation via the GI by 
numerical simulation including the dust self-gravity. Section |5] is devoted for conclusion and 
discussion. 



2. Equations and model 
2.1. Equations 

We consider a small region around the midplane which is rotating with the Keplerian 
frequency f2 at a distance r from a central star to study local dust motion and magnetohy- 
drodynamics. The coordinates that we use are {x, y, z) where x is the radial distance from 
r, y is tangential distance, and z is vertical distance from the disk midplane. 

We include centimeter to meter-size dust particles as super-particles. Total number of 
the super-particles is O(IO^) and each super-particle represents O(10^)-0(10'') small dust 
particles. The equation of motion of the i-th particle is given by 

= -2n X Vi + SQ^Xiic (v, - u) - V$, (1) 

dt Tf 

where u is the gas velocity at the location of the i-th particle, which is interpolated using 
gas velocities at the neighbor grid points. The third term in the r. h. s. (right hand side) 
represents the specific gas drag force to the dust particle. We adopt the Epstein drag force, 

Tf = ppa/{pgC,), (2) 

where Pp and a are the internal density and radius of dust particles, Pg and Cg are the spatial 
density and sound velocity of surrounding gas. The Epstein law is valid for centimeter to 
meter-size dust particles at ~ 5AU if gas colum n density is les s than twice as much as that 



of in the minimum solar nebula model (MMSN; lHayashilll98ll ). Even if we consider higher 
column density, the deviation in the drag force strength from Epstein law would not be 
significant. The last term in Equation ([T]) is the self-gravitational force of the dust particles. 
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We calculate the gravitational potential $ from the interpolated dust spatial density p^, 
solving the Poisson equation, 



= A-KGpd- 



(3) 



We calculate the self-gravity of the dust particles only in the situations where the GI is 
expected. We neglect the vertical gravity of the host star that ca uses settling of dust particles 
onto a thin layer (~ O.Olif where H is the disk scale height; ISchrapler &: Henningl 120041 ) 
near the midplane, because we are interested in radial concentration of dust but not vertical 
settling. While this would not affect the growth of MRI around the disk mid-plane that we 
simulate (|2;| < 0.25if), we do not have to resolve the thin dust layer, avoiding expensive 
computational cost. 

For the disk gas, we use the isothermal resistive MHD equations, 



dvi 

'dt 



9pg 
dt 



+ (u- V)u 



1 / 
V P+ — 

Pg V Stt 

—(3csQ:x. 



u-(v)). 



[B ■ V) B - 217 X u + 3fi2a;x 



+ V-(P3U) = 0, 



OB 

'dt 
P 



V X [(u X B) - r/(V X B)] 



Pgi 



(4) 

(5) 

(6) 
(7) 



where we assume constant Cg (an isothermal disk), (v) is the mean velocity field of the dust 
particles in the grid cell, e = Pdl Pg is the dust to gas ratio, and the term, — /3csf2x, expresses 
the global pressure gradient, which is separated from the local one. Let Pq and 5P be the 
global pressure and the local pressure variation due to the effect of MRI (P = Pq + ^P)- 
Assuming that Pq cx r"^, the global pressure gradient is given by 



}_dPo 

Pa dr 



^Pn 

Pg 



-q = qcg^l = — /3csf2, 

r 



(8) 



where H is the disk scale height defined hj H = Cg/Q. In our local model, (3 = qH/r 
is approximated to be constant. Note that both q and P are negative. Due to the radial 
pressure gradient, the gas rotation angular velocity is deviated from Keplerian one, as 



Q 1 



H 

r 



dlnP 
dlnr 




1 


'H' 


^ dlii6P\ 


2 


r 


dlnr I 



(9) 



The last term in the r. h. s. of Eq. is the dust drag force (back-reaction of gas drag 
force on the dust particles). Except for this term, the equations of motions for disk gas and 
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dust particles are the same as those in Paper II. The purpose of this paper is to study the 
effect of this term on dust concentration. When dust particles are accumulated and e takes 
a large value of > 0(1), this term would influence the gas flow. The induction equation 
(Eq. [B]) has the diffusion term by ohmic dissipation. Since we are interested in relatively 
inner regions (< a few dozens AU), we neglect the ambipolar diffusion that could influence 
MRI growth in the outer regions of > lOOAU (jChiang fc Murray-Clayl 120071 ). We treat the 
magnetic resistivity as a constant parameter for simplicity, though in reality it depends on 
the density of tiny grains (< fim). 



We scale length, time, and velocity by H, 1/fl, and c.^ in simulations. We solve 
the MHD equations by combining CIP scheme (lYabe fc Aokilll99ll ) and MOC-CT method 
(IStone &: NormanI Il992l ). The dust density is allocated to the closest eight grids in the 
three-dimensional space using cloud-in-cell (CIC) model. This algorithm st rictly conserves 
angula r momentum transfer from dust to gas by using a similar method as iJohansen et al. 
(120071 ). The boundary conditions in all directions are periodic. For the radial boundary, 
however, we take into account Keplerian differe ntial rotation with the shearing box model 
(jWisdom fc Tremaind Il988l : iHawley et al.l 119951 ) . While we are not interested in vertical 
sedimentation, we want to keep total dust mass in a whole simulation area, so that we 
adopt the periodic boundary condition also for vertical direction. The Poisson equation ([3]) 
is calculated by Fast Fourier Transform (FFT). This method requires the periodic bound- 
ary. For radial direction, according to the sheared boundary, we shift the phase azimuthally 
in Fourier spac e after the Fourier transform in the periodic azimuthal direction, following 
Johansen et al.l (120071 ). Our simulations are performed by a vector computer, NEC SX-9 at 
ISAS/JAXA. 



2.2. Initial conditions 

The setup for the simulation in this paper is the same as CASE2 in Paper II. We 
assume non-uniform to set marginal MRI state where localized dead (stable) and active 
(unstable) regions co-exist in the initial conditions. The initial magnetic field is Bq = 
(0, Bq sin 6, Bq cos 6), where 6 = 9{x) is the angle between the magnetic field and the vertical 
axis. We assume a constant value of Bq that is determined by plasma beta = 400, to establish 
the initial equilibration. With the constant magnetic resistivity of r? = 0.0 02ff we adopt. 



a threshold vertical magnetic field for MRI is -Bz^crit ~ 0.2i?o (jJinI Il996l ). We set radially 
varying 9 in which ^ = 0° [B^ ^ -Bz,crit) in the central zone and 9 = 85° [Bz < -Bz,crit) in the 
side regions (see Figure 2b in Paper II), such that MRI grows only in the central zone. We 
set the radial width of the initially active region as = lAH in all cases, and that of the 
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dead regions as = 4.0if in model-s40-* and Lg = 0.5H in model-s05-*, where the asterisk 
* represents other simulation parameters (see below). 

We assume equal-size dust particles with tj = 1.0/Q in model-s*-tlO-* except for tj = 
0.1 /Q in model-s40-t01-e010, neglecting their coalescence and fragmentation. At 5AU in 
MMSN, TfQ = 0.1 and 1.0 correspond to the dust sizes of approximately 3 and 30 centimeters, 
respectively. The super-particles are distributed uniformly such that the initial dust-to-gas 
ratio eo = 0.1 for model-s*-t*-e010 and eo = 0.01 for model-s40-tl0-e001. Note that the dust- 
to-gas ratio in our simulation box corresponds to that near the midplane layer. If vertically 
global sedimentation of dust particles is taken into account, the dust-to-gas ratio in our 
simulation box should be larger than that averaged over the whole disk. For comparison, we 
also present the results of Paper II without the dust drag as model-*-*-test. The self-gravity 
of dust particles in Equation ([T]) are switched on only in the saturated state in which the GI 
is expected. 

All of our simulations start with uniform gas density and pressure. The global pressure 
gradient is set to be /3 = —0.04. This assumed value of \(3\ is a few t imes smaller than that 



expected at ~ 5AU in MMSN. In order to compare the results with iJohansen et al.l (120071 ) 
and Paper II, we adopt the small value. As was argued in Paper II, the results would not be 
affected significantly by the value of The initial angular velocities of gas and particles 
are Uy/cg = — (3/2) (x/H) + (3/2 and Vy = — (3/2) (x/H), respectively. Because gas rotates 
slower, the particles migrate inward (negative direction of x) in the initial state. Initial 
disturbances are given randomly to the gas radial velocity with the amplitude of 0.001c<j. 
The size of our simulation box is {L^, Ly, L^) = ((2.5-9.5)iJ, 1.0H,0.5H) and the resolution 
is dx = dy = dz = O.OIH. Eight particles are distributed in each grid initially and the 
total number is ~ O(IO^). We have tested different number of distributed particles with 
different mass such that total mass is conserved and found that the results are converged if 
the distributed number of particles for each grid is > 8. The initial setup is summarized in 
Table E 



3. Effect of the dust drag on dust concentration 

In the simulations with dust drag, quasi-steady state is formed and the dust particles are 
concentrated around the outer-edge of the super-Keplerian region by inhomogeneous MRI. 
Here, we study the effect of the dust drag on the dust concentration by comparing the results 
with those without the dust drag. 
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3.1. TfQ, = 1.0 in a weak remnant turbulence 

In this subsection, we discuss the results with TfQ = 1.0, model-s40-tl0-e010 with 
eo = 0.10 and model-s40-tl0-e001 with eo = 0.01. Strong dust concentration was found 
in corresponding models without the dust drag in Paper II (model-s40-tl0-test), which is 
summarized in Figure [TJ MRI is excited only in the initially active central region (—0.71 < 
x/H < 0.71). The MRI turbulence transfers mass and angular momentum of disk gas to 
establish a local rigid rotation in the initially active zone through the turbulent viscosity. 
Then, in the outer half region of the active zone, gas flow is accelerated and migrates outward 
because of angular momentum gain. On the other hand, in the inner half region, it is 
decelerated and migrates inward. As a result, gas mass is moved from the central zone to 
the side zones and gas pressure is lowered in the central zone. The effect of the pressure 
modulation extends by radial scale of ~ H. Equation ([9]) shows that 

6u, = ^^^^ ^ ^ + = _o.02 + (10) 

This equation shows that super-Keplerian regions are associated with locally positive pressure 
gradient with some off-set due to global pressure gradient. We find that super-Keplerian 
regions are created in < a;/if < 2.0 and —2.8 ^ x/H < —2.5 (panel a and b). Although 
MRI turbulence has decayed in the result at tQ = 70, the magnetic field has not been 
dissipated in the central zone (panel c). The MRI is suppressed by the disappearance of 
shear motion. If the rigid rotation is perturbed toward the original Keplerian motion, the 
retrieved shear motion causes MRI again to recover the rigid rotation. Thus, this profile is 
stable and strong. 

Figure |2] presents the dust density in the saturated state in model-s40-tlO-e010 with 
Co = 0.10 (panel a) and model-s40-tl0-e001 with eo = 0.01 (panel b), in comparison with 
model-s40-tl0-test without the dust drag (panel c). Because the dust drag depends on spatial 
density of the dust particles, the results depend on eo- In all cases, after the particles are 
swept out of the active region by the MRI turbulence, they accumulate at the outer-edge of 
the super-Keplerian zone at x/iJ ^ 2.0 and -2.5. Particles leaving the simulation box from 
the small x (left hand) boundary reenter the simulation region from the large x (right hand) 
boundary after the shearing box correction is taken into account. The dust that reentered 
from the right hand boundary is halted at x/H ~ 2.0, resulting in further increasing of the 
dust density. The number of locations of dust concentration is fewer in the case with the dust 
drag, because the drag smoothes out small amplitude fluctuations of gas pressure. In the 
case with the drag, the individual dust concentrated areas are broader in model-s40-tl0-e010 
than in model-s40-tl0-e001, which is discussed below. We also found that velocity dispersion 
is lower for a denser clump, which was not observed in the case without the dust drag. The 
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effect of the reduced velocity dispersion will be discussed in section 4. 

Figure [3] shows the time evolution of the maximum density of dust particles (pd) in 
the simulation cells. The dust density is scaled by the gas density averaged in the whole 
simulation region ((pg)). The results are compared with those without the dust drag (dashed 
lines; model-Ls40-tlO-test) for eo = 0.10 and eo = 0.01. In the case without the dust drag, 
only concentration relative to the initial state is measured, so these lines are drawn by 
the evolution of concentration assuming eo = 0.10 or cq = 0.01. The maximum density 
continues to increase monotonically in this case. However, the growth of the maximum 
value is saturated in both cases with the dust drag. The saturation is faster and the increase 
rate is slightly smaller in model-s40-tl0-e010 than in model-s40-tl0-e001. 

In order to explain the broadening of the dust accumulated region, in addition to the 
three-dimensional simulations, we performed the two-dimensional {x-z) version of the model- 
s40-tl0-e010, in which the effect of the dust drag on the radial migration of the dust particles 
is more clearly shown. In Figure H^, we plot the difference of gas angular velocity from Kepler 
angular velocity, 6uy = {uy — fkep) /cs, near the sub- and super-Keplerian boundary in the 
two-dimensional simulation. The velocities are averaged over the vertical direction. At 
tVt = 55.0 (dashed line), the dust particles are expected to assemble at x/iJ ~ 1.75, where 
6uy = 0. At tQ = 87.4 (solid line), however, the region with 6uy ~ becomes broader 
(1.7 ^x/H < 1.9), because more dust particles have migrated to this region and their drag 
makes the gas flow close to Keplerian. Consequently, the dust particles are more broadly 
distributed there (Figure |Dd). Figure Ht schematically illustrates the effect of the dust drag. 
The radial velocity (oc 5uy) of a migrating dust particle becomes slower as the dust particle 
approaches the dust-concentrated region of 5uy = 0, like "traffic jam." Due to the modulation 
by the dust drag, the radial width of the Keplerian region is expanded, and the dust particles 
stop their inward migration before they reach the location at which 6uy = originally. Thus, 
the maximum dust density is self- regulated as shown in Figure O 

In the three-dimensional simulation, similar results are obtained, although small ampli- 
tude oscillations remain. The similarity implies that geometry is not the main cause for the 
broadening of the dust concentration region. 

3.2. Tf^l = 0.1 in a weak remnant turbulence 

In Paper II, for smaller particles with TfQ = 0.1 (model-s40-t01-test), we found that the 
dust concentration is not significant because the smaller dust particles are affected more by 
the turbulent diffusion. In the results in section 3.1, we found that the dust drag suppresses 
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the local turbulence and velocity dispersion of the dust particles. It could enhance the dust 
accumulation by deca ying the turbulence around the dust accumulated RrGcl, clS found in 



Johansen et al.l ( 120071 ) . 



However, the time evolution of the maximum scaled density of dust particles (p^/ (Pg)) in 
model-s40-t01-e010 is not significantly different from that in model-s40-t01-test (Figure [5]). 
In both cases, after pd rapidly increases by the diffusing-out from the initially active region 
by the MRI turbulence tVL ~ 25, it gradually increases and is saturated at ti7 > 50 — 70. The 
velocity dispersion of the dust particles is actually reduced from that in the case without the 
dust drag, but the effect is not strong enough to enhance the dust concentration. 



3.3. TfVt = 1.0 in a strong remnant turbulence 

In Paper II, we found that with the smaller initial dead region Lg = 0.55H, the viscosity 
in the saturated state is a ~ 10~^, which is much larger than a ~ 10~^ for the runs with 
Lg = 4.0H in model-Ls40-*. We found in Paper I that MRI turbulence does not decay if the 
magnetic Elssaser number radially averaged over the simulation region is smaller than unity 
in the initial state. Elsasser number is defined by 

Am,i,ve = vlz/V^, (11) 

where v\z = B^/ y^Anpg is z component of Alfven velocity and pg is spatial density of the 
disk gas. The run with Lg = 0.55H corresponds to Am ave ~ 0.5. 

The stronger remnant turbulence limited pd/ (pg) to the values less than 100 even for 
Tffl = 1.0 in the case without the dust drag. We performed the run with Ls = 0.55H, adding 
the dust drag (model-s05-tl0-e010). Although the drag force reduces the turbulent diffusion 
in the local concentrated region, the maximum dust density is not enhanced from that in 
the case without the dust drag (Figure E]). 



4. Planetesimal formation 

As shown in the previous section, the dust accumulation is self-regulated by the effect 
of the dust drag. The suppressed maximum density of dust particles is unfavorable for the 
gravitational instability (GI), while their reduced velocity dispersion is favorable. In this 
section, we analyze the results of previous runs without the self-gravity of dust particles to 
examine the possibility of subsequent GI. In the results of some favorable runs, we re-perform 
the simulations, including self-gravity among the dust particles, to demonstrate that the GI 



actually occurs. 



4.1. Analysis of gravitational instability 



The GI is expected to arise when the self-gravity of dust particles is dominant over their 
thermal fluctuation (velocity dispersion), in other words, when the radius (size) i? of a dust 
clump with mass M is smaller than the Jeans (Bondi) radius. 



:i2i 



where M = M{R) and cr = (t{R) are the total mass and mean velocity dispersion of particles 
in the region with distance R from the center of the clump, M^, is the host star's mass, Cg = 
VtH, and Vl = ^jGM^/r'^. The background shear is included in a. In order to numerically 
resolve GI, we set grid size in the x-direction such that dx < 0.5-Rj. 

The condition o f GI is often described by l inear analysi s of a uniform axisymmetric disk 



(e.g., 



Goldreich & Ward 


1973 


1 ( 


Toomrel 


1964 


), 



1 < Q 



TT 



(13) 



where Hd is the unperturbed solid column density. If we use M ~ TrS^i?^ and a ~ RQ, the 
Toomre's condition is identical to R < Rj except a numerical factor of 2. Because significant 
radial inhomogeneity develops before GI occurs in our case, we use the condition R < Rj 
that can be locally applied, rather than Q < 1. 

The mass of each super-particle m is given by Pdo/no, where pdo and Uq are the spatial 
density and number density of particles in the initial conditions. Then, the dust clump mass 
is given by 



Nr 

M = tuNr = — / 
no 



87r3 Nr 



H 

— I eo 

r 



f 1 M. 



(14) 



where Nr is the total particle number in R, Nro = (47r/3)no-R^ is its initial value, pdo = eoPgo 



and Pgo = Sgo/ v27iH by the assumption of a vertically isothermal disk. From Equtaions ( fT2l) 
and dH]), the scaled Jeans radius is given by 
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The initial dust-to-gas density ratio eo is given. The simulation results give dust enhancement 
in a clump Nfi/NjiQ, the clump size R/H, and velocity dispersion in the clump cr/cg (the 
equations of motions we use are normalized by H and Cg). From the results and simulation 
parameters, we can evaluate Rj for any given values of Hg^r"^ /M^, and H/r that are specified 
by disk model through Eq. (fT5|) . 

4.2. Effect of dust drag force on i?j 

The onset of GI depends on Nji and a as shown by Equation fllSI) . The dust drag 
suppresses Nji and reduces a (section [3]). Here, using the arguments in section 4.1., we 
examine the possibility of the GI in model-s40-tl0-e010, which is one of the most promising 
runs for the GI. 

Figure [7] shows Nr/Nrq from the densest grid point (panel a) and the corresponding 
a/cg (panel b) in the results of model-s40-tl0-e010 and model-s40-tl0-test at tVt = 58.0. 
Panel c shows Jeans radius -Rj calculated for M* = Mq, r = 5AU, H/r = 0.055 and 
Sgo = 150gcm~'^(~ Smmsn at r = 5AU, where Smmsn is gas column density in MMSN). 
While the particle concentration is lowered by the dust drag only slightly (panel a), the 
velocity dispersion is significantly reduced (panel b). Since the positive effect for the GI 
(reduction in the velocity dispersion) is dominated over the negative one (suppression in 
particle concentration), i?j calculated by Eq. (fT5|) is higher and the condition for the GI 
(i?j > R) is satisfied in the case with the dust drag (panel c). 

4.3. Simulation with dust self-gravity 

We carried out additional simulations including the self-gravity force of dust particles to 
demonstrate the formation of gravitationally bound clumps that may lead to planetesimals, 
in model-s40-tl0-e010. We set = Mq, r = 5AU, H/r = 0.055 and Sgo = 280gcm-3 ~ 
2Smmsn('" = 5AU). To reduce simulation cost, we introduced the self-gravity at tfl = 96 
when the dust concentration becomes saturated and Rj is much larger than the grid size dx. 

Figure [S] shows the time evolution of the gravitational collapse. Shortly after intro- 
duction of the self-gravity, the elongated high density region is kinked {tQ = 99) and it is 
separated into several clumps {tQ = 100). The clumps grow by accreting surrounding dust 
particles and other clumps [tQ = 120-140). 

To confirm that the clumps are gravitationally bound and estimate the mass of formed 
planetesimals, we define the range of a clump by its Hill's radius. The time evolution of the 
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clump is shown in Figure HI Figure shows the Hill's radius, where the grid size is repre- 
sented by a dashed line. Immediately after the introduction of the self-gravity at tQ = 96, 
the Hill radius exceeds the grid size. After that, the clump is numerically resolved. Figure [Hb 
shows that the velocity dispersion of the particles in the clump is always smaller than the 
surface escape velocity Vesc of the clump, which imphes that the clump is gravitationally 
bound. If the gravitational collapse continues, it may form a planetesimal, although this 
simulation does not have resolution to follow the subsequent collapse. 

Figure [nt shows the temporal development of the clump mass M, which may corre- 
spond to the planetesimal mass. The abrupt jumps at tfl ~ 102 and ~ 104 are caused by 
coUisional merging with other clumps. Since the destruction process is not properly included 
in our simulation, such rapid growth may be unrealistic. A conservative estimate for the 
planetesimal mass may be the mass before the abrupt jumps, that is, ~ 4 times Ceres mass. 
However, note that this mass is close to the resolution of our s imulation (Figure [9^) and the 



clump mass may be smaller in a higher- resolution simulation (jjohansen et al.ll2010l ). 



We also performed the simulation with the self-gravity in model-s40-tl0-e001, in which 
the Jeans radius is slightly larger than our grid size only for short interval. The GI is not 
found in this case, but it might be seen in a high-resolution simulation. 



4.4. Critical gas column density for gravitational instability 

In the simulation with addition of the self-gravity in section 4.3, we assumed Hgo that 
corresponds to ~ 2Smmsn at r = 5AU. On the other hand, simulations before adding 
self-gravity are scaled by 'Lg^r'^ /M^, and H/r. 

Here, fixing r = 5AU and H/r = 0.055, we apply the results of individual runs for 
various S^q to derive a sufficient condition for gas column density to cause the GI. The 
conditions for the GI is i? < Rj{R). Since Rj oc T,go (Equation [H]), the condition is more 
easily satisfied for larger S^o- Iii the saturated state in model-s40-tl0-e010, the condition is 
satisfied even at S^q/Smmsn ~ 1, while we showed the results with S^o = 2Smmsn in section 
4.3. For smaller eo (model-s40-tl0-e001), the critical column density is Sgo/SMMSN ~ 3. 
The smaller dust particles with Tffl = 0.1 have no chance to excite the GI even in the 
weak residual turbulence (model-s40-t01-e010) as long as Sc,o/Smmsn < 20. In the stronger 
remnant turbulence (Lg = 0.5H), the GI is not expected unless Sgo/SMMSN > 10 (model- 
s05-tl0-e010). 

Note that Rj is a function of eoS^o (Eq. [I5]). That is, the possibility of the GI depends 
on the total column density of dust particles, but not on eo. Thus, for example, i?j should 



- 15 - 



be similar between the result with Sgo/SMMSN = 1 in model- Ls40-tl0-e010 (eo = 0.10) and 
Sgo/^MMSN = 10 in model-Ls40-tlO-e001 (eo = 0.01) at the same r. From the results of 
simulations in this paper, it is inferred that the GI may occur when p^ccrit ^ 0.03pg^MMSN- 

5. Conclusion and Discussion 

We have studied the dust concentration including the " dust drag force" onto gas (back- 
reaction of the gas drag exerted onto the dust particles) in a quasi-steady state created by 
inhomogeneous MRI found by Paper I and II, by performing the three-dimensional resistive 
MHD simulation including dust particles as super-particles. Since the inertia of the particles 
is taken into account, the dust drag force modulates gas flow in the dust concentrated regions. 
We also examined the possibility of the planctcsimal formation via gravitational instability 
(Gl) with analysis using Jeans radius of dense dust regions and performed simulations with 
adding self-gravity of the dust particles to demonstrate that gravitationally bound clumps 
are actually formed by the GI. 

If MRI active and dead zones initially coexist, mass and angular momentum transfer 
associated by non-uniformly growing MRI turbulence changes the slightly sub-Keplerian gas 
flow in the initial state to a quasi-steady MRI-stable state in which super- and sub-Keplerian 
regions are radially adjacent to each other (Paper I), and the dust particles are concentrated 
at the outer edge of the super-Keplerian region (Paper II). In this paper, we found that 
the introduction of the dust drag broadens the dust accumulated regions while it reduces 
velocity dispersion of the particles, depending on the turbulent level and the friction time of 
the dust particles. We found that the positive effect (the reduction in velocity dispersion) is 
generally dominated over the negative effect (the broadening of the dust accumulated region). 
Consequently, in the case with dust drag, the Gl is expected in the case of weak remnant 
turbulence (the turbulent viscosity a ~ 10~^) and meter-size particles with r/Q = 1.0, if the 
initial dust spatial density is a few times larger than that of MMSN. The GI is regulated by 
the absolute value of the dust spatial density, but not by the dust-to-gas ratio. 

Smaller dust particles {rfVt = 0.1) are also less likely to cause the Gl even in the 
weak remnant turbulence case, because they are more strongly coupled with gas turbulent 
motion. Since dust particles should have size distribution, the spatial density contributed 
from meter-size particles must be a few times larger than total dust density of MMSN for the 

GI. However, dust settling increases dust density in the layer of midplane that corresponds to 
our simulation box and it could compensate for the effect of size distribution. If the vertical 
dust distribution is Gaussian (oc e^^^/^^^) before the dust settling and it is assumed that 
most of dust particles settle down to our simulation box with Lz = 0.5H, the averaged dust- 
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to-gas ratio of our simulation box is enhanced by a factor of five from the initial dust-to-gas 
ratio of the whole disk. 

In the models for dust growth in turbulent eddies proposed by other authors, high colli- 
sion velocity between the dust particles excited by the turbulence may result in fragmentation 
rather than coalescence, which is not favored for planetesimal formation. However, in our 
model, MRI turbulence is almost terminated after it transforms the gas flow to the quasi- 
steady state with the pressure bump, so that the collision velocities between dust particles 
are as small as < 0.6-0. 7m/s at r = 3-5AU, which can avoid fragmentation at mutual 
collisions. 

Radially non-uniform excitation of MRI is an essential point for the emergence of the 
pressure bump in our model. The growth rate of MRI depends on the strength of the 
magnetic field and resistivity. In Paper II, we found non-uniform resistivity produces the 
same quasi-steady state as non-uniform magnetic field that we assume in this paper. In 
section 1, we raised a possibility of formation of the active zone radially sandwiched by 
dead zones due to non-uniform resistivity near the snow line. However, the location of 
the snow line and dead zones are coupled with disk evolution due to viscous diffusion and 
photoevaporation and also with growth, fragmentation and migration of dust particles. Thus, 
to evaluate the possibility of radially "local" planetesimal formation proposed by this paper, 
full-scale coupled evolution of the snow line, the dead zone, the disk gas ionization degree 
and dust growth needs to be studied theoretically and by observation with ALMA. 

We thank for detailed comments by an anonymous referee. This work was supported 
by Grant-in- Aid for JSPS Fellows (208778). The simulations presented in this paper were 
performed by NEC SX-6 at ISAS/JAXA. 
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Fig. 1. — Results of model-s40-tl0-test described in Paper II. Time evolution of vertically 
averaged values of (a) pressure P, (b) gas angular velocity Uy and (c) number density of 
particles n. P and n are normalized by the initial values [Pq and no), and Uy is normalized 
by sound speed c^. The dotted, dashed and bold lines represent the snapshots at tfl = 0,40 
and 70, respectively. The two vertical dotted-lines are the boundaries between the initially 
active (unstable) and dead (stable) regions. MRI is initially excited only in the region 
between the two dotted lines. 
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Fig. 2. — The dust density at the saturated state in (a) model-Ls40-tlO-e010 {rffl = 1.0 and 
eo = 0.10), (b) model-Ls40-tlO-e001 [TfQ = 1.0 and eo = 0.01) and (c) model-Ls40-tlO-test 
{TfQ = 1.0 without dust drag force). The samphng time is tfl = 104. The initially active 
region is located between the two white lines. 
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Fig. 3. — Time evolution of the maximum dust concentration in model-Ls40-tlO-e010, -eOOl 
and -test. The sohd and dashed hues represent the results of eo = 0.10 and eo = 0.01, 
respectively. All lines represent the dust density in the cell having the highest density in 
the whole simulation region, which is normalized by the gas density averaged over the whole 
region ((pg)). The thin dotted lines show the result without the dust drag (model-Ls40-tlO- 
test). In this result, only concentration relative to the initial state is measured, so these lines 
are drawn by assuming eo = 0.10 or eo = 0.01. 
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Fig. 4. — Broadening of the dust concentrated region by the dust drag, (a) The difference 
between the gas angular velocity and Kepler angular velocity and (b) the vertically averaged 
dust density at tQ = 55.4 (dashed lines) and tQ = 70.0 (solid lines) in model-Ls40-tlO- 
eOlO {TfQ = 1.0 and eo = 0.10). These figures are magnification of the area around the 
concentrated region where Uy = fkep- (c) Schematic illustration of "traffic jam" of the dust 
particles. 
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Fig. 5. — Same as Figure [3] but for model-s40-t01-e010 (bold solid line; TfQ = 0.1 and 
eo = 0.10) and -test (thin dotted line; Tffl = 0.1 without the dust drag force). 
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Fig. 6. — Same as Figure [3] but for model-s05-tl0-e010 (bold solid line; Lg = 0.55H) and 
-test (thin dotted line; without the dust drag force). 
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Fig. 7. — Estimation of possibility of the GI by Eq. ( IT5l) . In all panels, the solid and dashed 
lines represent model-s40-tl0-e010 and model-s40-tl0-test at tQ = 58.0, respectively, (a) 
The number Npt of particles within distance R from a densest grid normalized by the initial 
value iV/jo- (b) The velocity dispersion of the particles in R. (c) The radius Rj calculated 
for each R for M* = Mq, r = 5AU and S^o = 150gcm^^ ~ Smmsn(^ = 5AU). In the region 
over the dotted line {Rj > R), the GI is expected. 
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Fig. 8. — Simulation of the GI with self-gravity of the dust particles in model-Ls40-tlO-e010. 
The top panel shows the dust density at tQ = 96.0, at which the self-gravity of particles 
is added. The bottom panels show the time evolution of the GI. The different gray colors 
represent the isosurface of the dust-density \og{pd/ {pg)) = 1.0,2.0 and 3.0. The several 
clumps become bounded gravitationally. 
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Fig. 9. — Evolution of the gravitationally bounded clump in model-Ls40-tlO-e010. (a) The 
Hill's radius, (b) the surface escape velocity {vesc] the solid line) and the velocity dispersion of 
particles {a; the dashed line) in the Hill's radius, and (c) the total mass of the dust particles 
in the Hill's radius. 
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Run 


Ls 




eo 


self- 


-gravity 


model-s40-tl0-e010 


4.0// 


1.0 


0.10 


off 


and on 


-tlO-eOOl 


4.0// 


1.0 


0.01 


off 


and on 


-tlO-test 


4.0// 


1.0 


0.0 




off 


-tOl-eOlO 


4.0// 


0.1 


0.10 




off 


-tOl-test 


4.0// 


0.1 


0.0 




off 


model-s05-tl0-e010 


0.55// 


1.0 


0.10 




off 


-tlO-test 


0.55// 


1.0 


0.0 




off 



Tabic 1: Setup of individual runs. Lg is the radial width of initially dead region; r/ is friction 
time; eo is initial dust-to-gas density ratio. Model-s40-tl0-e010 and eOOl are also re-started 
with introduction of self-gravity of particles. 



